Spin and charge order in doped Hubbard model: long- wavelength collective modes 
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Determining the ground state properties of the two-dimensional Hubbard model has remained an outstanding 
problem. Applying recent advances in constrained path auxiliary-field quantum Monte Carlo techniques and 
simulating large rectangular periodic lattices, we calculate the long-range spin and charge correlations in the 
ground state as a function of doping. At intermediate interaction strengths, an incommensurate spin density 
wave (SDW) state is found, with antiferromagnetic order and essentially homogeneous charge correlation. The 
wavelength of the collective mode decreases with doping, as does its magnitude. The SDW order vanishes 
beyond a critical doping. As the interaction is increased, the holes go from a wave-like to a particle-like state, 
and charge ordering develops which eventually evolves into stripe-like states. 
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The Hubbard model is one of the simplest and most 
fundamental models for quantum many-body systems. Since 
the discovery of high-Tc superconductors in 1986, the two- 
dimensional (2D) repulsive Hubbard model has received great 
interest, as a potential model for describing the essential 
physics of the copper-oxygen plane in cuprates iQl- Thanks to 
intensive analytic and numerical investigations, some aspects 
of its phase diagram have been understood [3|]. For example, 
at half-filling (one electron per lattice site), the system has 
long-range antifeiTomagnetic (AF) order 10. H S Bl- How- 
ever, many basic issues still remain unknown or controversial. 
For instance, what happens to the AF order when the system 
is doped? This is important not only for understanding the 
magnetic properties, which the Hubbard model was originally 
designed to describe, but also in the context of high-Tc, which 
shares the same parameter regime and is beUeved to be closely 
related to AF fluctuations. 

The difficulty in treating the Hubbard model underscores a 
more general theme, namely the challenge of accurate treat- 
ment of strongly correlated systems. Recent experimental 
work ijsil with cold fermionic atoms in optical lattices offers 
a promising new avenue — potentially direct simulations of 
Hubbard-like models with "lattice emulators" |9]. We believe 
this increases, rather than decreases, the demand on "tradi- 
tional" numerical simulations. Although all numerical meth- 
ods have their limitations, high quality data on the Hubbard 
model will provide guidance and allow direct comparison with 
experiments, thereby creating a new level of synergy to tackle 
the problem of strong electron correlations. 

There are many earlier investigations of the spin and charge 
correlations in the 2D Hubbard model S H S S 0, E [ll . 
H El [S El E [B [iS El , 

but reaching large system sizes 
with sufficient accuracy has been a challenge. Recently, we 
have calculated the equation of state from low to intermediate 
interaction strengths on periodic lattices of up to 16 x 16 
|2C1. The constrained path Monte Carlo (CPMC) method 
1 2]J, I22II was generalized to incorporate a boundary condition 



integration technique 1231], which removed short-range finite- 
size effects. It was found that, immediately upon doping, the 
thermodynamic stability condition is violated. This implied 
the existence of a spatially inhomogeneous phase. In the ab- 



sence of long-range collective modes, the results are an accu- 
rate representation of the thermodynamic limit, and the insta- 
bility would indicate phase separation. On the other hand, 
if long-range collective modes existed whose characteristic 
length exceeds the size of the super cell ('^ 16), they would 
not be fully captured in the simulations. The nature of the AF 
fluctuation in the doped Hubbard model thus remains to be 
resolved. 

Here we address the question, by employing recent algo- 
rithmic advances |23, 23, 25] in CPMC and simulating rect- 



angular supercells on parallel computers. Much larger linear 
dimension (128) is reached than in previous studies {^^ 16), 
and detailed measurements are obtained of the spin-spin and 
charge-charge correlations in the ground state, at intermedi- 
ate interaction strengths where our method is very accurate. 
Our results show that long wavelength collective modes of in- 
commensurate spin density wave (SDW) states appear as the 
system is doped, with AF spin order but essentially homo- 
geneous charge correlation. Charge correlation develops as 
the interaction is further increased. We quantify the nature of 
such states, and discuss how they relate to the "stripe" states 
at large interactions. 

The Hamiltonian for the single-band Hubbard model is 

where ^ (cj.a) creates (annihilates) an electron with spin 
(7 ((7 =1,1) at lattice site j, and 6 connects all possible 
nearest-neighbor sites. The supercell has N ^ Lx x Ly sites. 
The density is n = [N^ +Ni)/N, where A^o- is the number 
of electrons with spin C7; doping is li = 1 — n. We imple- 



ment twist-averaged boundary conditions II23I1 . under which 
the wave function gains a phase when electrons hop around 



lattice boundaries: 



•) = 



(0v,0y) 



are ran- 



where L is the unit vector along L, and 
dom twist angles over which we average. 

The generalized CPMC method lHaEHIllIll] used here 
obtains the many-body ground state by repeated projections 
with e^'^^ (t is the projection time step), as in standard quan- 
tum Monte Carlo (QMC). The two-body part, e"'^^^ is de- 
coupled via the Hubbard-Stratonovich transformation into a 
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sum over one-body projectors in Ising fields 0]. The projec- 
tion is then reaUzed efficiently by importance-sampled ran- 
dom walks with non-orthogonal Slater determinants (SDs), 
where the one-body projectors propagate one SD into another, 
and the many-dimensional sum over Ising fields is performed 
by Monte Carlo. The usual fermion sign/phase problem is 
controlled approximately by a global phase condition on the 
SDs. This is the only approximation in the method. The ba- 
sic idea of the approximation is as follows. The many-body 
ground state is given by |^o) = where \(p) are 

SDs sampled by the QMC, and their probability distribution 
will give w(0) (> 0). Because the Schrodinger equation is 
linear, is degenerate with — jl'o)- A. trivial effect in a 
deterministic representation, this can cause the determinants 
{(j)), in a random walk, to move back and forth between the 
two sets of solutions. In a simulation, precisely when a |0) 
turns from one to the other can not be detected, because the 
continuous stochastic evolution of the orbitals can lead to an 
exchange without any two orbitals ever overlapping. This is 
the sign problem. We use a trial wave function \^t) to make 
the detection, by requiring {^t\^) > 0. Because each |0) is a 
full many-electron wave function, the sign of its overlap with a 
\^t) is expected to be quite insensitive to the details of \^t)- 

In extensive benchmarks in Hubbard models lEol 21, 3] 
as well as in atoms, molecules, and solids ^J^, this general 
framework has demonstrated accuracy equaling or surpassing 
the most accurate (non-exponential scaling) many-body com- 
putational methods available. In the Hubbard model, the en- 
ergy is typically within < 0.5% of the exact diagonalization 
results for U = 4t ll20ll . At half-filling where the approxima- 
tion is the most severe (with free-electron l^r)), the method 
gives an energy per site of —0.8559(4) for the infinite lat- 
tice, compared to the estimated exact result of —0.8618(2) 
1 15l 12811 . (The method can be made exact at half-filling by 
removing the constraint lElll .) 

In order to probe correlations at long range, we study 
rectangular supercells of 4 x Ly, 8 x Ly and 16 x L,,. The 
largest system size simulated in this work is N ^ 1024, us- 
ing ^^(lOOO) processors on the Cray XT4 supercomputer The 
first indication of a long wavelength collective mode is seen 
in the ground-state energy. Figure [T] shows how the energy 
per site, e, varies as Ly is increased. Each energy has been 
averaged over 20-1000 0-values, and all controllable QMC 
biases (e.g.. Trotter and population size 12111 ') have been re- 
moved. The error bars are estimated by combining statistical 
error and 0-point fluctuation. Twist-averaging eliminates ki- 
netic energy finite-size effects (shell and lattice size) 1I20I] . At 
/i = 1/4 for example, we have reconfirmed that the energy re- 
mained essentially constant as Ly was varied from 8 to 64. 

At h = 1/16, the energies remain above the line from 
Maxwell construction up to Ly ~ 16, where it shows a signifi- 
cant drop and falls below the line. In the inset, the equation of 
state (EOS), e(n) vs. n, is shown for 8x8. The EOS is con- 
cave for n e («e, 1) ||2^. The critical density ric is determined 
by the Maxwell construction, which gives a phase separation 
line tangent to the EOS: eM(«) ~ h/ hc£{nc) + {1 ~ h/ hc)e{n — 
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FIG. 1 : Energy variations with linear dimension. The QMC ground 
state energy per site, e (in units of /), is shown for U = 4t for a 
sequence of supercell sizes 8 x Ly. Blue squares are with an FE \^t}, 
while red circles are with a UHF Ifj)- The main graph is for h = 
1/16. The green line is determined by Maxwell construction, which 
is illustrated in the inset. To magnify the vertical scale in the inset, a 
linear shift of £M(n) has been applied, so the tangent line is horizontal 
(dashed orange), and the EOS is plotted as e(«) — eM(")- 



1). The drop for Ly > 16 indicates that the instability occurs 
only in smaller supercells, in which a state with long-range 
correlation is frustrated. 

We use two different types of \^t), to help gauge the effect 
of the constraint. The first is the free-electron (FE) wave func- 
tion, which is of course homogeneous with no long-range cor- 
relation. The second is the unrestricted Hartree-Fock (UHF) 
solution, which has broken spatial symmetry and static long- 
range spin and charge order We have verified that, in the 
paramagnetic phase below /i^, the two types lead to statisti- 
cally indistinguishable QMC results. As seen in Fig. [T] the 
same is true at /j = 1 /16, except for large systems {Ly > 16), 
where the UHF \^t) gives lower energy. This is consistent 
with UHF being a better wave function in a system where an 
SDW can develop, as further discussed below. 

We calculate the spin-spin correlation function: 

'^^•W = 7?L(('W,T ~"r+r',i) («r',T -nr',i)) , (2) 
r' 

which measures the correlation between two spins separated 
by a lattice vector r = (/,,,/, ). The corresponding structure 
factor is 5,(k) = e''""C, (r). Similarly, we calculate the 
charge-charge correlation Q (r), defined by replacing "— " in 
Eq. (|2]i with "+", and its structure factor 5t: (k). 

The QMC results for G(r) are shown in Fig. |2] An AF 
correlation is seen clearly in the density plots, in which the 
signs alternate for near neighbors. The staggered correlation 
function: C',{lx,ly) = (-l)'iCs(/i,/,,), is also plotted with sta- 
tistical error bars. The curves fall into two groups, for even 
and odd Ix, respectively. With perfect AF correlation, each 
group would be a constant function of ly. In these systems the 
AF correlation patterns are modulated by a wave along ly. A 
K phase-shift occurs at the nodes where C[{lx,ly) crosses zero. 
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FIG. 2: Spin-spin correlation function in 8 x 32 (top panel) and 
8 X 64 (bottom panel) lattices at = 1/16 and U = '^t. The results are 
averaged over 20 ©-points. UHF |'P7-)'s are used. The upper part in 
each panel is a 3D density plot (color theme in the upper right corner 
of the graph). The lower part of each panel shows the staggered 
correlation, with curves of different colors representing different (f 's. 
Due to symmetry, only G [0,Lv/2] is shown. The range of Cj(r) is 
restricted to [—0.1,0.1], so the "self-peak" near the origin is cut off. 



The wave essentially doubles as L, is doubled, with compara- 
ble SDW amplitude. We see that the wavelength of the spin 
modulation is A 32 at this density, consistent with Fig. [T] 
where the energy lowering only occurs at L, ^16. 

We next study the spin-spin correlation as a function of dop- 
ing. Calculations are done at three densities, h = 3/32, 1/16, 
and 1 /32, respectively, for a 4 x 64 lattice, at U ^ At. The 
staggered spin-spin correlation function ^^(r) are shown in 
Fig. [3] where the modulation and n phase shifts are clearly 
seen. The wavelength of the modulation decreases with dop- 
ing, beginning at half-filling {h = 0) where the wavelength is 
X ~oa. (Quantitatively, our data is consistent with the wave- 
length being inversely proportional to h, although statistical 
error bars on X are large.) The strength of Cj(r) in the incom- 
mensurate SDW state also decreases with doping, and appears 
to vanish at a critical value of /ze ~ 0. 15 ± 0.05, where the sys- 
tem turns into a paramagnetic liquid. 

The results in Fig.|2]were obtained using UHF I'I't'), while 
those in Fig.|3]were from FE |^r)- The consistency between 
them is reassuring. To generate the UHF |^r), we used the 
minimum U for which a UHF solution exists ({/ ~ 1.3-1.5f 
for /i = 1/16 and ^ 3.5f for /i = 1 /4), in order to minimize 
the effect of broken translational invariance. A weak static 
long-range order is present in the UHF solution. For example, 
in the 8 x 64 system, calculated from the UHF \^t) itself 
was ^(5 x 10^'*). We see that this was enhanced by a factor 
of 200 in the QMC, to ^(0.1). On the other hand, at large 
distance the variation in the UHF charge-charge correlation 
Cc was ^(1.5 X 10""^), which remained in the QMC 
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FIG. 3: Doping dependence of the long wavelength incommensurate 
SDW state. The staggered spin-spin correlation function C'^{lxJy) 
is plotted vs. /,., at three different densities. Calculations are done 
using free-electron (FE) \^t}- The system is a 4 x 64 supercell, with 
U = 4t and &/7t = (-0.8410, -0.9198). Colors label different value 
of /v's as in Fig.|2] As doping is increased, the wavelength of the 
modulating wave decreases, as does the amplitude of the SDW. 



(zero within error bars). With the FE \^t), the SDW struc- 
ture in Fig.^emerged spontaneously. QMC results with UHF 
\^t) always showed long-range order, while we do find vari- 
ations when the FE \^t) is used: for some 0's no long-range 
SDW order is seen, only "short"-range incommensurate AF 
correlations. In such cases, the calculated QMC energy tends 
to be higher than that using UHF \^t)- In contrast, when a 
long-range SDW is also seen with the FE \^t), the QMC en- 
ergies are more consistent between the two types of |^7')'s. 
We interpret this as the system favoring long-range order. The 
rectangular supercells break the symmetry between the x- and 
y- directions. In the thermodynamic limit a combination of 
the linear SDWs may be present. 

The spin correlation discussed above {U = At) is not ac- 
companied by charge inhomogeneities. The effect of stronger 
interactions is examined in Fig. |4] which displays results for 
f//f = 4, 8 and 12, with doping of /? = 1 /16 in a 4 x 32 lat- 
tice. The spin structure factor 5v(k) is plotted along the line 
cut k = {K,ky). At = 4f a pronounced peak can already 
be seen at {n, 15;r/16), consistent with the spin-spin correla- 
tion in Fig. |2] (The split of the ^^(k) peak with doping had 
also been observed in earlier simulations Jsl S [3 111 ) As 
U is increased, the peak value increases rapidly. The charge 
structure factor is plotted along (0,^:,,). Except for the trivial 
peak at the origin, ^^ (k) is broad with little features alU = 4f . 
Above {/ = 8f, a peak appears at (0,;r/8), indicating the de- 
velopment of a charge-charge correlation. 

In the inset of Fig. HI the real-space density profile is plotted 
along (0, ly). AtU ~ At, results using FE and UHF (generated 
with U = I. At) trial wave functions both give a constant den- 
sity. At larger U, the same UHF j^r) (from U — I. At) turned 
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FIG. 4: Spin-spin and charge-charge structure factors as a func- 
tion of interaction strength U. Results are shown for a 4 x 32 su- 
percell, with @/k = (-0.8410,-0.9198) and at doping h = 1/16. 
FE and UHF indicate the type of I'I't-) used in the QMC calculation. 
The inset in the right panel shows real space density profile with 
the same color coding as in the main figures. S.s(k) has a peak at 
k = {7t, 15;r/16), with the peak value growing rapidly with increas- 
ing U. 5e(k), which has no feature atU = 4t, shows a corresponding 
peak at k = (0, ;r/8) for large U. The real-space density is a constant 
atU = 4t, but develops periodic modulation at larger U. 

out to be sufficient to "pin" the many-body solution into a 
broken translational symmetry charge density wave state. The 
density profile provides a way to visualize the nature of the 
state. At [/ = I2t, the region of maximum density tends to sat- 
urate at p = 1, while "stripes" appear at the boundaries which 
separate AF spin domains with a n phase shift. This is con- 
sistent with density matrix renormalization group results lfl6ll 
of stripe states in Hubbard ladders for [/ > 8f . The real-space 
characteristic length of the charge correlation is 1/2 that of 
the spin correlation, as Fig. |4] shows. These results suggest 
that, at intermediate U, holes are in a uniform "liquid" state 
with no long-range correlation, while at the large U limit they 
enter a "Wigner-crystal" state forming stripes. 

The difficulty in treating the Hubbard model arises from 
the multiple competing energy scales separated by tiny differ- 
ences. The system can fall into one phase or the other due 
to a small bias in the calculation or in simulation boundary 
condition. The challenge for numerical calculations is to min- 
imize the effect of such biases (intrinsic accuracy, system size, 
etc). Our calculations reach much larger systems than possi- 
ble otherwise. In this work, we have carefully removed biases 
other than the effect of the constrained path approximation. 
To address the latter, we have used trial wave functions with 
opposite properties (uniform FE vs. broken-symmetry UHF) 
to examine the robustness and consistency of the results. 

To conclude, we have presented numerical results from 
constrained path QMC to characterize the magnetic proper- 
ties in doped 2D Hubbard model. At intermediate interaction 
strengths U /t ^ A, the ground state has incommensurate anti- 
ferromagnetic SDW order with long wavelength modulation. 



The wavelength of the SDW and the strength of the spin or- 
der both decrease with doping, and the state vanishes below a 
critical density, when the system enters a paramagnetic "liq- 
uid" phase. In the SDW state there is essentially no charge 
correlation, with the holes in a wave-like state. As U in- 
creases, accompanying charge correlation develops, with the 
holes becoming localized at the nodal positions of the modu- 
lating wave. Thus in the strong interaction regime {U >~ lOf) 
the system evolves into a stripe-like state. Many topics remain 
for future work, including quantitative aspects of these states 
and their implications on superconductivity. 
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support by ARO (56693-PH) and NSF (DMR-0535592). 
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